 
module iondrag
  !-------------------------------------------------------------------------------
  ! Purpose:
  !   Calculate ion drag tendency and apply to horizontal velocities.
  !   Also calculate joule heating tendency and apply to neutral temperature. 
  ! 
  ! Subroutines:
  !   iondrag_init (initialize module)
  !   iondrag_calc (calculate ion drag tensors)
  !   iondrag_tend (ion drag tendency)   
  !   qjoule_tend (joule heating)
  !
  ! Calling sequence:
  !   inti
  !     iondrag_init
  !   tphysac
  !     iondrag_calc
  !       iondrag_tend
  !       qjoule_tend
  !
  ! Dependencies:
  !   Magnetic field from apex module
  !   ExB ion drifts from exbdrift module
  !
  ! Author:
  !   B. Foster Feb, 2004
  !
  !-------------------------------------------------------------------------------

  use shr_kind_mod ,only: r8 => shr_kind_r8
  use ppgrid       ,only: pcols, pver, begchunk, endchunk
  use cam_history  ,only: addfld, add_default, outfld, horiz_only
  use physics_types,only: physics_state, physics_ptend, physics_ptend_init
  
  use physics_buffer, only : pbuf_get_index, physics_buffer_desc, pbuf_get_field
  use perf_mod     ,only: t_startf, t_stopf
  use cam_logfile  ,only: iulog

  use interpolate_data, only: lininterp
  use spmd_utils,       only: masterproc

  use phys_control,  only: waccmx_is
  use cam_abortutils,only: endrun
  use time_manager,  only: is_first_step

  implicit none

  save

  private                         ! Make default type private to the module

  !-------------------------------------------------------------------------------
  ! Public interfaces:
  !-------------------------------------------------------------------------------
  public :: iondrag_register         ! Register variables in pbuf physics buffer
  public :: iondrag_init             ! Initialization
  public :: iondrag_calc             ! ion drag tensors lxx,lyy,lxy,lyx
  public :: iondrag_readnl
  public :: iondrag_timestep_init
  public :: iondrag_inidat
  public :: do_waccm_ions

  interface iondrag_calc
     module procedure iondrag_calc_ions
     module procedure iondrag_calc_ghg
  end interface

  !-------------------------------------------------------------------------------
  ! Private data
  !-------------------------------------------------------------------------------

  ! Namelist variables
  character(len=256) :: efield_lflux_file = 'coeff_lflux.dat'
  character(len=256) :: efield_hflux_file = 'coeff_hflux.dat'
  real(r8)           :: efield_potential_max = huge(1._r8) ! max cross cap potential kV
  logical            :: empirical_ion_velocities = .true.

  real(r8),parameter :: amu    = 1.6605387e-27_r8  ! atomic mass unit (kg)

  integer  :: ntop_lev = 1
  integer  :: nbot_lev = 0
  integer  :: id_xo2, id_xo1             ! indices to tn and major sp
  integer  :: id_o2p, id_op, id_nop      ! indices to ions
  integer  :: id_elec, id_n

  !Physics buffer indices
  integer :: PedConduct_idx  = 0
  integer :: HallConduct_idx = 0
  integer :: ui_idx      = 0 ! index to zonal drift from edynamo
  integer :: vi_idx      = 0 ! index to meridional drift from edynamo
  integer :: wi_idx      = 0 ! index to vertical drift from edynamo
  integer :: indxTe  = -1    ! pbuf index for electron temperature
  integer :: indxTi  = -1    ! pbuf index for ion temperature

  logical  :: xo2_slvd, xo1_slvd, o2p_slvd, op_slvd, nop_slvd

  real(r8) :: rmass_op                   ! mass of O+ (g/mole)
  real(r8) :: rmass_o2p                  ! mass of O2+ (g/mole)
  real(r8) :: rmass_nop                  ! mass of NO+ (g/mole)
  real(r8) :: rmass_o1                   ! mass of O (g/mole)
  real(r8) :: rmass_o2                   ! mass of O2 (g/mole)
  real(r8) :: rmass_n2                   ! mass of N2 (g/mole)

  !-------------------------------------------------------------------------------
  ! Inverted masses (for multiply in loops rather than divide):
  !-------------------------------------------------------------------------------
  real(r8) :: rmi_o1
  real(r8) :: rmi_o2
  real(r8) :: rmi_n2
  real(r8) :: rmi_op
  real(r8) :: rmi_o2p
  real(r8) :: rmi_nop
  real(r8) :: rmi_op_kg
  real(r8) :: rmi_o2p_kg
  real(r8) :: rmi_nop_kg

  ! GHG
  !-------------------------------------------------------------------------

  ! Private data
  integer, parameter :: plevtiod = 97   

  real(r8) alamxx(plevtiod)
  real(r8) alamxy(plevtiod)

  real(r8) pshtiod(plevtiod)           ! TIME pressure scale height
  real(r8) pshccm(pver)                ! CCM pressure scale height

  real(r8) alamxxi(pver)              ! alamxx interpolated to waccm grid
  real(r8) alamxyi(pver)              ! alamxy interpoalted to waccm grid

  logical doiodrg
  logical, protected :: do_waccm_ions = .false.

  !
  ! Data statement for ALAMXX
  !
  data alamxx /                                                     &
       0.13902E-17_r8, 0.22222E-17_r8, 0.34700E-17_r8, 0.53680E-17_r8, 0.83647E-17_r8, &
       0.13035E-16_r8, 0.20254E-16_r8, 0.31415E-16_r8, 0.48944E-16_r8, 0.75871E-16_r8, &
       0.11584E-15_r8, 0.17389E-15_r8, 0.25786E-15_r8, 0.37994E-15_r8, 0.58088E-15_r8, &
       0.95179E-15_r8, 0.19052E-14_r8, 0.47869E-14_r8, 0.14284E-13_r8, 0.45584E-13_r8, &
       0.14756E-12_r8, 0.48154E-12_r8, 0.14844E-11_r8, 0.39209E-11_r8, 0.83886E-11_r8, &
       0.14213E-10_r8, 0.20304E-10_r8, 0.27449E-10_r8, 0.39276E-10_r8, 0.59044E-10_r8, &
       0.83683E-10_r8, 0.11377E-09_r8, 0.14655E-09_r8, 0.19059E-09_r8, 0.28338E-09_r8, &
       0.46326E-09_r8, 0.73966E-09_r8, 0.11785E-08_r8, 0.18789E-08_r8, 0.31037E-08_r8, &
       0.53919E-08_r8, 0.97251E-08_r8, 0.17868E-07_r8, 0.33041E-07_r8, 0.61265E-07_r8, &
       0.11406E-06_r8, 0.20912E-06_r8, 0.39426E-06_r8, 0.76691E-06_r8, 0.15113E-05_r8, &
       0.29545E-05_r8, 0.55644E-05_r8, 0.97208E-05_r8, 0.16733E-04_r8, 0.28101E-04_r8, &
       0.36946E-04_r8, 0.44277E-04_r8, 0.50982E-04_r8, 0.57526E-04_r8, 0.64190E-04_r8, &
       0.71471E-04_r8, 0.80311E-04_r8, 0.96121E-04_r8, 0.11356E-03_r8, 0.14131E-03_r8, &
       0.18695E-03_r8, 0.26058E-03_r8, 0.36900E-03_r8, 0.50812E-03_r8, 0.66171E-03_r8, &
       0.80763E-03_r8, 0.92583E-03_r8, 0.10038E-02_r8, 0.10382E-02_r8, 0.10333E-02_r8, &
       0.99732E-03_r8, 0.93994E-03_r8, 0.86984E-03_r8, 0.79384E-03_r8, 0.71691E-03_r8, &
       0.64237E-03_r8, 0.57224E-03_r8, 0.50761E-03_r8, 0.44894E-03_r8, 0.39624E-03_r8, &
       0.34929E-03_r8, 0.30767E-03_r8, 0.27089E-03_r8, 0.23845E-03_r8, 0.20985E-03_r8, &
       0.18462E-03_r8, 0.16233E-03_r8, 0.14260E-03_r8, 0.12510E-03_r8, 0.10955E-03_r8, &
       0.95699E-04_r8, 0.83347E-04_r8/

  !
  ! Data statement for ALAMXY
  !
  data alamxy /                                                     &
       0.74471E-24_r8, 0.22662E-23_r8, 0.69004E-23_r8, 0.20345E-22_r8, 0.58465E-22_r8, &
       0.16542E-21_r8, 0.46240E-21_r8, 0.12795E-20_r8, 0.35226E-20_r8, 0.96664E-20_r8, &
       0.26650E-19_r8, 0.76791E-19_r8, 0.25710E-18_r8, 0.10897E-17_r8, 0.56593E-17_r8, &
       0.30990E-16_r8, 0.16792E-15_r8, 0.85438E-15_r8, 0.40830E-14_r8, 0.18350E-13_r8, &
       0.79062E-13_r8, 0.33578E-12_r8, 0.13348E-11_r8, 0.45311E-11_r8, 0.12443E-10_r8, &
       0.27052E-10_r8, 0.49598E-10_r8, 0.86072E-10_r8, 0.15807E-09_r8, 0.30480E-09_r8, &
       0.55333E-09_r8, 0.96125E-09_r8, 0.15757E-08_r8, 0.25896E-08_r8, 0.48209E-08_r8, &
       0.96504E-08_r8, 0.18494E-07_r8, 0.34296E-07_r8, 0.61112E-07_r8, 0.10738E-06_r8, &
       0.18747E-06_r8, 0.32054E-06_r8, 0.52872E-06_r8, 0.83634E-06_r8, 0.12723E-05_r8, &
       0.18748E-05_r8, 0.26362E-05_r8, 0.36986E-05_r8, 0.52079E-05_r8, 0.72579E-05_r8, &
       0.98614E-05_r8, 0.12775E-04_r8, 0.15295E-04_r8, 0.18072E-04_r8, 0.20959E-04_r8, &
       0.19208E-04_r8, 0.16285E-04_r8, 0.13628E-04_r8, 0.11784E-04_r8, 0.11085E-04_r8, &
       0.11916E-04_r8, 0.14771E-04_r8, 0.20471E-04_r8, 0.29426E-04_r8, 0.42992E-04_r8, &
       0.62609E-04_r8, 0.90224E-04_r8, 0.12870E-03_r8, 0.18281E-03_r8, 0.26029E-03_r8, &
       0.37224E-03_r8, 0.53254E-03_r8, 0.75697E-03_r8, 0.10623E-02_r8, 0.14660E-02_r8, &
       0.19856E-02_r8, 0.26393E-02_r8, 0.34473E-02_r8, 0.44327E-02_r8, 0.56254E-02_r8, &
       0.70672E-02_r8, 0.88174E-02_r8, 0.10960E-01_r8, 0.13613E-01_r8, 0.16934E-01_r8, &
       0.21137E-01_r8, 0.26501E-01_r8, 0.33388E-01_r8, 0.42263E-01_r8, 0.53716E-01_r8, &
       0.68491E-01_r8, 0.87521E-01_r8, 0.11196E+00_r8, 0.14320E+00_r8, 0.18295E+00_r8, &
       0.23321E+00_r8, 0.29631E+00_r8/

  logical :: ionvels_read_from_file = .false.

contains

!==============================================================================     

  subroutine iondrag_register
!-----------------------------------------------------------------------
! Register E and B fields.
!
! Register iondrag variables with physics buffer:
!
! Hall and Pedersen conductivities
! 
! pcols dimension and lchnk assumed here
!
!-----------------------------------------------------------------------
    use exbdrift,        only: exbdrift_register
    use physics_buffer,  only: pbuf_add_field, dtype_r8

    ! E and B fields
    call exbdrift_register()

    if ( waccmx_is("ionosphere") ) then
       ! Pedersen Conductivity and Hall Conductivity
       call pbuf_add_field('PedConduct',  'physpkg', dtype_r8, (/pcols,pver/), PedConduct_idx )
       call pbuf_add_field('HallConduct', 'physpkg', dtype_r8, (/pcols,pver/), HallConduct_idx)
    end if

  end subroutine iondrag_register

!================================================================================================

  subroutine iondrag_readnl(nlfile)

    use namelist_utils,  only: find_group_name
    use units,           only: getunit, freeunit
    use spmd_utils, only: mpicom, masterprocid, mpicom, mpi_character, mpi_logical, mpi_real8

    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input

    ! Local variables
    integer :: unitn, ierr
    character(len=*), parameter :: subname = 'iondrag_readnl'

    namelist /iondrag_nl/ efield_lflux_file, efield_hflux_file, efield_potential_max, empirical_ion_velocities

    if (masterproc) then
       unitn = getunit()
       open( unitn, file=trim(nlfile), status='old' )
       call find_group_name(unitn, 'iondrag_nl', status=ierr)
       if (ierr == 0) then
          read(unitn, iondrag_nl, iostat=ierr)
          if (ierr /= 0) then
             call endrun(subname // ':: ERROR reading namelist')
          end if
       end if
       close(unitn)
       call freeunit(unitn)

    end if

    call mpi_bcast (efield_lflux_file, len(efield_lflux_file), mpi_character, masterprocid, mpicom, ierr)
    call mpi_bcast (efield_hflux_file, len(efield_hflux_file), mpi_character, masterprocid, mpicom, ierr)
    call mpi_bcast (efield_potential_max, 1, mpi_real8, masterprocid, mpicom, ierr)
    call mpi_bcast (empirical_ion_velocities, 1, mpi_logical,  masterprocid, mpicom, ierr)

    if (masterproc) then
       write(iulog,*) 'iondrag options:'
       write(iulog,*) '  efield_lflux_file = '//trim(efield_lflux_file)
       write(iulog,*) '  efield_hflux_file = '//trim(efield_hflux_file)
       write(iulog,*) '  efield_potential_max = ',efield_potential_max
       write(iulog,*) '  empirical_ion_velocities = ',empirical_ion_velocities
    end if

  end subroutine iondrag_readnl

  !================================================================================================

  subroutine iondrag_init( pref_mid )
    use constituents, only: cnst_get_ind
    use short_lived_species, only: slvd_index

    !-------------------------------------------------------------------------------
    ! Iondrag initialization, called from inti.F90.
    !-------------------------------------------------------------------------------


    !-------------------------------------------------------------------------------
    ! dummy arguments
    !-------------------------------------------------------------------------------
    real(r8), intent(in) :: pref_mid(pver)

    integer :: k, err

    integer :: cnst_ids(7)

    doiodrg = .false.
    do_waccm_ions = .false.

    !-------------------------------------------------------------------------------
    ! find lower bnd for iondrag
    !-------------------------------------------------------------------------------
    if( pref_mid(1) < 0.1_r8 ) then
       do k = 1, pver
          if (pref_mid(k) < 50._r8) nbot_lev = k
       end do
    end if

    if (nbot_lev > 0) then
       doiodrg = .true.
    endif

    if ( .not. doiodrg .and. masterproc ) then
       write(iulog,*) ' '
       write(iulog,*) 'iondrag_init: Does not have waccm level. Ion drag does not apply. '
       write(iulog,*) ' '
       return
    endif

    if( masterproc ) then
       write(iulog,*) ' '
       write(iulog,*) 'iondrag_init: nbot_lev,press = ',nbot_lev,pref_mid(nbot_lev)
       write(iulog,*) ' '
    end if

    call cnst_get_ind( 'e', id_elec, abort=.false. )
    if (id_elec < 0) then
       id_elec = slvd_index( 'e' )
    endif
    call cnst_get_ind( 'Op', id_op, abort=.false. )
    if (id_op < 0) then
       id_op = slvd_index( 'Op' )
       if (id_op > 0) then
          op_slvd = .true.
       endif
    else
       op_slvd = .false.
    endif
    call cnst_get_ind( 'O2p', id_o2p, abort=.false. )
    if (id_o2p < 0) then
       id_o2p = slvd_index( 'O2p' )
       if (id_o2p > 0) then
          o2p_slvd = .true.
       endif
    else
       o2p_slvd = .false.
    endif
    call cnst_get_ind( 'NOp', id_nop, abort=.false. )
    if (id_nop < 0) then
       id_nop = slvd_index( 'NOp' )
       if (id_nop > 0) then
          nop_slvd = .true.
       endif
    else
       nop_slvd = .false.
    endif
    call cnst_get_ind( 'O', id_xo1, abort=.false. )
    if (id_xo1 < 0) then
       id_xo1 = slvd_index( 'O' )
       if (id_xo1 > 0) then
          xo1_slvd = .true.
       endif
    else
       xo1_slvd = .false.
    endif
    call cnst_get_ind( 'O2', id_xo2, abort=.false. )
    if (id_xo2 < 0) then
       id_xo2 = slvd_index( 'O2' )
       if (id_xo2 > 0) then
          xo2_slvd = .true.
       endif
    else
       xo2_slvd = .false.
    endif
    call cnst_get_ind( 'N', id_n, abort=.false. )
    if (id_n < 0) then
       id_n = slvd_index( 'N' )
    endif
 
    cnst_ids = (/ id_elec, id_op, id_o2p, id_nop, id_xo1, id_xo2, id_n /)

    if ( all( cnst_ids > 0 ) ) then
       do_waccm_ions = .true.
    endif

    if ( do_waccm_ions ) then
       call ions_init
    else
       call ghg_init(pref_mid)
    endif

    if (.not.doiodrg) return

    ! Add to masterfield list
    call addfld('UIONTEND',(/ 'lev' /),'A','M/S2','u-tendency due to ion drag')
    call addfld('VIONTEND',(/ 'lev' /),'A','M/S2','v-tendency due to ion drag')

!
! Indices to 3d ion drifts :
    ui_idx = pbuf_get_index('UI')
    vi_idx = pbuf_get_index('VI')
    wi_idx = pbuf_get_index('WI')

    indxTe = pbuf_get_index( 'TElec',errcode=err )
    indxTi = pbuf_get_index( 'TIon',errcode=err )

  end subroutine iondrag_init

  !================================================================================================
  subroutine ions_init

    use efield,       only: efield_init
    use exbdrift,     only: exbdrift_init
    use mo_chem_utls, only: get_spc_ndx
    use chem_mods,    only: adv_mass
    use mo_chem_utls, only: get_inv_ndx
    use chem_mods,    only: fix_mass
    use phys_control, only: phys_getopts

    !-------------------------------------------------------------------------------
    ! local variables
    !-------------------------------------------------------------------------------
    integer  :: id

    logical :: history_waccm

    call phys_getopts(history_waccm_out=history_waccm)

    !-------------------------------------------------------------------------------
    ! initialize related packages: electric field
    !-------------------------------------------------------------------------------

    call efield_init (efield_lflux_file, efield_hflux_file, efield_potential_max)
    call exbdrift_init( empirical_ion_velocities )

    id = get_spc_ndx('Op')
    rmass_op  = adv_mass(id)
    id = get_spc_ndx('O2p')
    rmass_o2p = adv_mass(id)
    id = get_spc_ndx('NOp')
    rmass_nop = adv_mass(id)
    id = get_spc_ndx('O')
    rmass_o1  = adv_mass(id)
    id = get_spc_ndx('O2')
    rmass_o2 = adv_mass(id)
    id = get_inv_ndx('N2')
    rmass_n2 = fix_mass(id)

    rmi_o1     = 1._r8/rmass_o1
    rmi_o2     = 1._r8/rmass_o2
    rmi_n2     = 1._r8/rmass_n2
    rmi_op     = 1._r8/rmass_op
    rmi_o2p    = 1._r8/rmass_o2p
    rmi_nop    = 1._r8/rmass_nop
    rmi_op_kg  = 1._r8/(rmass_op *amu)
    rmi_o2p_kg = 1._r8/(rmass_o2p*amu)  
    rmi_nop_kg = 1._r8/(rmass_nop*amu)

    !-------------------------------------------------------------------------------
    ! Set up fields to history files.
    !-------------------------------------------------------------------------------
    call addfld('BNORTH',  horiz_only,  'I', 'nT', 'northward component of magnetic field (nT)')
    call addfld('BEAST' ,  horiz_only,  'I', 'nT', 'eastward component of magnetic field (nT)')
    call addfld('BDOWN' ,  horiz_only,  'I', 'nT', 'downward component of magnetic field (nT)')
    call addfld('BMAG'  ,  horiz_only,  'I', 'nT', 'magnetic field magnitude (nT)' )

    call addfld('QIONSUM',   (/ 'lev' /), 'I','S-1' ,'Ion prod sum')
    call addfld('ELECDEN',   (/ 'lev' /), 'I','CM-3','NE (ion sum)')
    call addfld('SIGMAPED',  (/ 'lev' /), 'I', 'siemens/m', 'Pederson conductivity' )
    call addfld('SIGMAHAL',  (/ 'lev' /), 'I', 'siemens/m', 'Hall conductivity' )
    call addfld('LAMDA1'  ,  (/ 'lev' /), 'I' ,'S-1','LAMDA PED')
    call addfld('LAMDA2'  ,  (/ 'lev' /), 'I' ,'S-1','LAMDA HALL')

    call addfld('LXX',       (/ 'lev' /), 'I','S-1','LXX')
    call addfld('LYY',       (/ 'lev' /), 'I','S-1','LYY')
    call addfld('LXY',       (/ 'lev' /), 'I','S-1','LXY')
    call addfld('LYX',       (/ 'lev' /), 'I','S-1','LYX')

    !
    ! Joule heating, and tn before and after joule heating tendencies are applied:
    !
    call addfld( 'QJOULE',   (/ 'lev' /), 'I', 'K/s' , 'Joule Heat' )  ! joule heating
    if (history_waccm) then
       call add_default( 'QJOULE  ', 1, ' ' )                                ! joule heating (K/s)
    end if
!
! 3d drifts from either edynamo or exbdrift.
!
    if (empirical_ion_velocities) then
      call addfld('UI',(/ 'lev' /),'I','m/s', 'UI Zonal empirical ExB drift from exbdrift')
      call addfld('VI',(/ 'lev' /),'I','m/s', 'VI Meridional empirical ExB drift from exbdrift')
      call addfld('WI',(/ 'lev' /),'I','m/s', 'WI Vertical empirical ExB drift from exbdrift')
    endif

  end subroutine ions_init

  !========================================================================

  subroutine ghg_init (pref_mid)

    !
    ! initialization for ion drag calculation
    !

    !------------------Input arguments---------------------------------------

    real(r8), intent(in) :: pref_mid(pver)           ! model ref pressure at midpoint   

    !-----------------local workspace---------------------------------------
    integer k
    integer kinv

    real(r8) rpsh                                    ! ref pressure scale height

    real(r8), parameter :: preftgcm = 5.e-5_r8       ! TIME GCM reference pressure (Pa)

    !------------------------------------------------------------------------

    ! With the defualt values of nbot_lev and ntop_lev, ion drag calcualtion are NOT carried out
    nbot_lev=0
    ntop_lev=1

    do k = 1, pver
       rpsh=log(1e5_r8/pref_mid(k))
       if (rpsh .gt. 14._r8) nbot_lev  = k
    end do
    if (nbot_lev .gt. ntop_lev) doiodrg=.true.
    if (masterproc) then
       write(iulog,fmt='(a15)') 'From IONDRAGI:'
       write(iulog,fmt='(1a12,1i10)') 'NTOP_LEV  =',ntop_lev
       write(iulog,fmt='(1a12,1i10)') 'NBOT_LEV  =',nbot_lev
       write(iulog,*) 'IONDRAG flag is',doiodrg
    endif
    if (.not.doiodrg) return

    !     obtain TIME/GCM pressure scale height
    pshtiod(1)=-17._r8
    do k=2,plevtiod
       pshtiod(k)=pshtiod(k-1)+0.25_r8
    enddo

    !     map TIME-psh into CCM-psh 
    pshtiod=pshtiod-log(preftgcm/1E5_r8)

    !     CCM psh
    !     note that vertical indexing is inverted with respect to CCM standard
    do k=1,pver
       kinv=pver-k+1
       pshccm(kinv)=log(1e5_r8/pref_mid(k))
    enddo

    !     vertical interpolation 
    write(iulog,*) ' '
    write(iulog,*) 'iondragi: before lininterp for alamxx'
    write(iulog,*) '          nlatin,nlatout =',plevtiod,pver
    write(iulog,*) '          yin'
    write(iulog,'(1p,5g15.8)') pshtiod
    write(iulog,*) '          yout'
    write(iulog,'(1p,5g15.8)') pshccm
    write(iulog,*) ' '

    call lininterp (alamxx  ,pshtiod,plevtiod, alamxxi   ,pshccm,pver)

    call lininterp (alamxy  ,pshtiod,plevtiod, alamxyi   ,pshccm,pver)

    !     invert indeces back to CCM convention
    alamxxi(1:pver)=alamxxi(pver:1:-1)
    alamxyi(1:pver)=alamxyi(pver:1:-1)

    return

  end subroutine ghg_init

  !================================================================================================

  subroutine iondrag_timestep_init
    use efield, only: get_efield

    if (do_waccm_ions) then
      if (empirical_ion_velocities .or. (is_first_step().and..not.ionvels_read_from_file)) then
       ! Compute the electric field
       call t_startf ('efield')
       call get_efield
       call t_stopf ('efield')
      endif
    endif

  end subroutine iondrag_timestep_init

  !================================================================================================
  subroutine iondrag_calc_ions( lchnk, ncol, state, ptend, pbuf,  delt )
    !-------------------------------------------------------------------------------
    ! Calculate ion drag tensors lxx,lyy,lxy,lyx.
    ! Also calculate Pedersen and Hall conductivities.
    ! This is called from tphysac.
    !-------------------------------------------------------------------------------

    use mo_apex,only: & ! (pcols,begchunk:endchunk)
         bnorth,      & ! northward component of magnetic field (nT)
         beast,       & ! eastward component of magnetic field (nT)
         bdown,       & ! downward component of magnetic field (nT)
         bmag           ! magnetic field magnitude (nT)
    use physconst,     only: avogad, boltz
    use chemistry,     only: imozart
    use mo_mean_mass,  only: set_mean_mass
    use exbdrift,      only: exbdrift_ion_vels
    use short_lived_species, only: slvd_pbf_ndx => pbf_idx
    use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_set_field

    !-------------------------------------------------------------------------------
    ! dummy arguments
    !-------------------------------------------------------------------------------
    integer,intent(in)   :: lchnk               ! current chunk index
    integer,intent(in)   :: ncol                ! number of atmospheric columns
    real(r8), intent(in) :: delt                ! time step (s)
    type(physics_state), intent(in), target    :: state ! Physics state variables
    type(physics_ptend), intent(out)   :: ptend ! Physics tendencies
    
    type(physics_buffer_desc), pointer :: pbuf(:)

    !-------------------------------------------------------------------------------
    ! Local variables
    !-------------------------------------------------------------------------------
    integer :: i,k ! loop indices

    real(r8) :: sqrt_te                              ! sqrt(te)
    real(r8) :: sqrt_tnti                            ! sqrt(tnti)
    real(r8) :: wrk
    real(r8),parameter :: dipmin = 0.17_r8           ! minimum dip angle (tuneable)
    real(r8),parameter :: emass  = 9.1093819e-31_r8  ! electron mass (kg)
    real(r8),parameter :: qe     = 1.6021765e-19_r8  ! electronic charge (coulombs)
    real(r8),parameter :: colfac = 1.5_r8            ! collision factor (tuneable)
    real(r8),parameter :: boltzmann  = 1.e7_r8 * boltz
    real(r8),parameter :: avo    = avogad*1.e-3_r8 ! (molecules/mole)
    ! real(r8),parameter :: rmass_op  = 15.9989_r8     ! mass of O+
    ! real(r8),parameter :: rmass_o2p = 31.9983_r8     ! mass of O2+
    ! real(r8),parameter :: rmass_nop = 30.0056_r8     ! mass of NO+
    ! real(r8),parameter :: rmass_o1  = 16._r8         ! mass of O
    ! real(r8),parameter :: rmass_o2  = 32._r8         ! mass of O2
    ! real(r8),parameter :: rmass_n2  = 28._r8         ! mass of N2

    !-------------------------------------------------------------------------------
    ! Inverted masses (for multiply in loops rather than divide):
    !-------------------------------------------------------------------------------
    ! real(r8),parameter :: rmi_o1     = 1._r8/rmass_o1
    ! real(r8),parameter :: rmi_o2     = 1._r8/rmass_o2
    ! real(r8),parameter :: rmi_n2     = 1._r8/rmass_n2
    ! real(r8),parameter :: rmi_op     = 1._r8/rmass_op
    ! real(r8),parameter :: rmi_o2p    = 1._r8/rmass_o2p
    ! real(r8),parameter :: rmi_nop    = 1._r8/rmass_nop
    ! real(r8),parameter :: rmi_op_kg  = 1._r8/(rmass_op *amu)
    ! real(r8),parameter :: rmi_o2p_kg = 1._r8/(rmass_o2p*amu)  
    ! real(r8),parameter :: rmi_nop_kg = 1._r8/(rmass_nop*amu)

    real(r8), target :: tn(pcols,pver) ! neutral gas temperature (deg K)
    real(r8) :: xo2     (pcols,pver)   ! O2  (mmr)
    real(r8) :: xo1     (pcols,pver)   ! O   (mmr)
    real(r8) :: xn2     (pcols,pver)   ! N2  (mmr)
    real(r8) :: o2p     (pcols,pver)   ! O2+ (mmr)
    real(r8) :: op      (pcols,pver)   ! O+  (mmr)
    real(r8) :: nop     (pcols,pver)   ! NO+ (mmr)
    real(r8) :: barm    (pcols,pver)   ! mean molecular weight (g/mole)
    real(r8) :: xnmbar  (pcols,pver)   ! for unit conversion to volume density
    real(r8) :: tnti    (pcols,pver)   ! average of tn and ti
    real(r8) :: o2_cm3  (pcols,pver)   ! o2 volume density (cm-3)
    real(r8) :: o1_cm3  (pcols,pver)   ! o  volume density (cm-3)
    real(r8) :: n2_cm3  (pcols,pver)   ! n2 volume density (cm-3)
    real(r8) :: o2p_cm3 (pcols,pver)   ! O2+ (cm-3)
    real(r8) :: op_cm3  (pcols,pver)   ! O+  (cm-3)
    real(r8) :: nop_cm3 (pcols,pver)   ! NO+ (cm-3)
    real(r8) :: ne_sigmas(pcols,pver)  ! electron density for conductivities
    real(r8) :: lamda1  (pcols,pver)   ! sigped*b**2/rho
    real(r8) :: lamda2  (pcols,pver)   ! sighal*b**2/rho
    real(r8) :: lxxnorot(pcols,pver)   ! XX before rotation
    real(r8) :: lyynorot(pcols,pver)   ! YY before rotation
    real(r8) :: lxynorot(pcols,pver)   ! XY before rotation

    !-------------------------------------------------------------------------------
    ! Ion-neutral momentum transfer collision frequencies.
    ! rnu_xxx_xx = ratio of collision to gyro-frequences for O2+, O+, NO+.
    !-------------------------------------------------------------------------------
    real(r8) :: rnu_o2p_o2(pcols,pver)   ! O2+ ~ O2 collision freq (resonant, T dependent)
    real(r8) :: rnu_op_o2 (pcols,pver)   ! O+  ~ O2 collision freq (non-resonant)
    real(r8) :: rnu_nop_o2(pcols,pver)   ! NO+ ~ O2 collision freq (non-resonant)
    real(r8) :: rnu_o2p_o (pcols,pver)   ! O2+ ~ O  collision freq (non-resonant)
    real(r8) :: rnu_op_o  (pcols,pver)   ! O+  ~ O  collision freq (resonant, T dependent)
    real(r8) :: rnu_nop_o (pcols,pver)   ! NO+ ~ O  collision freq (non-resonant)
    real(r8) :: rnu_o2p_n2(pcols,pver)   ! O2+ ~ N2 collision freq (non-resonant)
    real(r8) :: rnu_op_n2 (pcols,pver)   ! O+  ~ N2 collision freq (non-resonant)
    real(r8) :: rnu_nop_n2(pcols,pver)   ! NO+ ~ N2 collision freq (non-resonant)
    real(r8) :: rnu_o2p   (pcols,pver)   ! [[o2p~o2]n(o2)+[o2p~o]n(o)+[o2p~n2]n(n2)]/w(o2p)
    real(r8) :: rnu_op    (pcols,pver)   ! [[op ~o2]n(o2)+[op ~o]n(o)+[op ~n2]n(n2)]/w(op )
    real(r8) :: rnu_nop   (pcols,pver)   ! [[nop~o2]n(o2)+[nop~o]n(o)+[nop~n2]n(n2)]/w(nop)
    real(r8) :: rnu_ne    (pcols,pver)   ! electron ~ neutral collision frequency (s-1)

    real(r8) :: press        (pcols)     ! pressure at interface levels (dyne/cm^2)
    real(r8) :: qe_fac       (pcols)     ! unit conversion factor for conductivities
    real(r8) :: dipmag       (pcols)     ! magnetic dip angle 
    real(r8) :: decmag       (pcols)     ! magnetic declination
    real(r8) :: btesla       (pcols)     ! magnetic field (teslas)
    real(r8) :: sindip       (pcols)     ! sin(dipmag)
    real(r8) :: sin2dip      (pcols)     ! sindip^2
    real(r8) :: sindec       (pcols)     ! sin(decmag)
    real(r8) :: cosdec       (pcols)     ! cos(decmag)
    real(r8) :: sin2dec      (pcols)     ! sindec^2
    real(r8) :: cos2dec      (pcols)     ! cosdec^2
    real(r8) :: omega_o2p    (pcols)     ! angular gyrofrequency for o2+ (s-1)
    real(r8) :: omega_op     (pcols)     ! angular gyrofrequency for o+ (s-1)
    real(r8) :: omega_nop    (pcols)     ! angular gyrofrequency for no+ (s-1)
    real(r8) :: omega_e      (pcols)     ! electron angular gyrofrequency (s-1)
    real(r8) :: omega_o2p_inv(pcols)     ! inverse of o2+ gyrofrequency
    real(r8) :: omega_op_inv (pcols)     ! inverse of o+  gyrofrequency
    real(r8) :: omega_nop_inv(pcols)     ! inverse of no+ gyrofrequency
    real(r8) :: omega_e_inv  (pcols)     ! inverse of electron gyrofrequency

    !-------------------------------------------------------------------------------
    ! Ion drag coefficients output:
    !-------------------------------------------------------------------------------
    real(r8) :: lxx(pcols,pver)  ! lambda XX coefficients (s-1)
    real(r8) :: lyy(pcols,pver)  ! lambda YY coefficients (s-1)
    real(r8) :: lxy(pcols,pver)  ! lambda XY coefficients (s-1)
    real(r8) :: lyx(pcols,pver)  ! lambda YX coefficients (s-1)

    !-------------------------------------------------------------------------------
    ! Conductivities output:
    !-------------------------------------------------------------------------------
    real(r8) :: sigma_ped (pcols,pver)   ! pedersen conductivity (siemens/m)
    real(r8) :: sigma_hall(pcols,pver)   ! hall conductivity (siemens/m)

    real(r8) :: qout(pcols,pver)         ! temp for outfld

    real(r8), dimension(:,:), pointer :: q_xo1,  q_xo2, q_o2p, q_op, q_nop

    real(r8), dimension(:,:), pointer :: tE  	 ! electron temperature in pbuf (K) 
    real(r8), dimension(:,:), pointer :: tI   	 ! ion temperature in pbuf (K) 

    if (.not.doiodrg) return

    call outfld ('BNORTH', bnorth(:,lchnk), pcols, lchnk )
    call outfld ('BEAST' , beast(:,lchnk),  pcols, lchnk )
    call outfld ('BDOWN' , bdown(:,lchnk),  pcols, lchnk )
    call outfld ('BMAG'  , bmag(:,lchnk),   pcols, lchnk )

    if ( xo1_slvd ) then
      call pbuf_get_field(pbuf, slvd_pbf_ndx, q_xo1, start=(/1,1,id_xo1/), kount=(/pcols,pver,1/) )
    else
      q_xo1 => state%q(:,:,id_xo1)
    endif
    if ( xo2_slvd ) then
      call pbuf_get_field(pbuf, slvd_pbf_ndx, q_xo2, start=(/1,1,id_xo2/), kount=(/pcols,pver,1/) )
    else
      q_xo2 => state%q(:,:,id_xo2)
    endif
    if ( o2p_slvd ) then
      call pbuf_get_field(pbuf, slvd_pbf_ndx, q_o2p, start=(/1,1,id_o2p/), kount=(/pcols,pver,1/) )
    else
      q_o2p => state%q(:,:,id_o2p)
    endif
    if ( op_slvd ) then
      call pbuf_get_field(pbuf, slvd_pbf_ndx, q_op, start=(/1,1,id_op/), kount=(/pcols,pver,1/) )
    else
      q_op => state%q(:,:,id_op)
    endif
    if ( nop_slvd ) then
      call pbuf_get_field(pbuf, slvd_pbf_ndx, q_nop, start=(/1,1,id_nop/), kount=(/pcols,pver,1/) )
    else
      q_nop => state%q(:,:,id_nop)
    endif

    !-------------------------------------------------------------------------------
    ! Define local tn and major species from state (mmr): 
    !-------------------------------------------------------------------------------
    do k = 1,pver
       do i = 1,ncol
          tn (i,k) = state%t(i,k)
          xo2(i,k) = q_xo2(i,k)              ! o2  (mmr)
          xo1(i,k) = q_xo1(i,k)              ! o   (mmr)
          xn2(i,k) = 1._r8 - (xo2(i,k) + xo1(i,k)) ! n2  (mmr)
          xn2(i,k) = max( 1.e-20_r8,xn2(i,k) )
          o2p(i,k) = q_o2p(i,k)              ! o2+ (mmr)
          op (i,k) = q_op(i,k)               ! o+  (mmr)
          nop(i,k) = q_nop(i,k)              ! no+ (mmr)
       end do
    end do

    !--------------------------------------------------------------------------------------------------
    ! For WACCM-X, grab electron (TE) and ion (TI) temperatures from pbuf from ionosphere module
    !--------------------------------------------------------------------------------------------------
    if ( indxTe>0 .and. indxTi>0 ) then
      call pbuf_get_field(pbuf, indxTe, tE)
      call pbuf_get_field(pbuf, indxTi, tI)
    else
       tE => tn
       tI => tn
    endif

    qout(:ncol,:) = o2p(:ncol,:) + op(:ncol,:) + nop(:ncol,:)
    call outfld ('QIONSUM ', qout, pcols, lchnk)

    !-------------------------------------------------------------------------------
    ! calculate empirical ExB drift velocities if not using edynamo drifts
    ! (if edynamo calculated the drifts, then they are already in pbuf%ui, etc.)
    !-------------------------------------------------------------------------------
    if (empirical_ion_velocities .or. (is_first_step().and..not.ionvels_read_from_file)) then
      call t_startf ( 'exbdrift_ion_vels' )
      call exbdrift_ion_vels( lchnk, ncol, pbuf)
      call t_stopf  ( 'exbdrift_ion_vels' )
    endif

    !-------------------------------------------------------------------------------

    do i = 1,ncol
       btesla(i) = bmag(i,lchnk)*1.e-9_r8 ! nT to teslas (see bmag in apex module)
       !-------------------------------------------------------------------------------
       ! Angular gyrofrequency of O+, O2+ and NO+ (s-1):
       !-------------------------------------------------------------------------------
       omega_op (i) = qe*btesla(i)*rmi_op_kg
       omega_o2p(i) = qe*btesla(i)*rmi_o2p_kg
       omega_nop(i) = qe*btesla(i)*rmi_nop_kg
       !-------------------------------------------------------------------------------
       ! Electron angular gyrofrequency (s-1):
       !-------------------------------------------------------------------------------
       omega_e(i) = qe*btesla(i)/emass 
       !-------------------------------------------------------------------------------
       ! Invert now, so we can multiply rather than divide in loops below:
       !-------------------------------------------------------------------------------
       omega_op_inv (i) = 1._r8/omega_op(i)
       omega_o2p_inv(i) = 1._r8/omega_o2p(i)
       omega_nop_inv(i) = 1._r8/omega_nop(i)
       omega_e_inv(i)   = 1._r8/omega_e(i)
       !-------------------------------------------------------------------------------
       ! Magnetic field geometry (used below in rotation of lambdas):
       !-------------------------------------------------------------------------------
       dipmag(i) = atan( bdown(i,lchnk)/sqrt(bnorth(i,lchnk)**2+beast(i,lchnk)**2) )
       decmag(i) = -atan2( beast(i,lchnk),bnorth(i,lchnk) )
       cosdec(i) = cos( decmag(i) )
       sindec(i) = sin( decmag(i) )
       if( abs(dipmag(i)) >= dipmin ) then
          sindip(i) = sin(dipmag(i))
       else
          if( dipmag(i) >= 0._r8 ) then
             sindip(i) = sin( dipmin )
          else
             sindip(i) = sin( -dipmin )
          end if
       end if
       sin2dip(i) = sindip(i)**2
       sin2dec(i) = sindec(i)**2
       cos2dec(i) = cosdec(i)**2
    end do

    ! write(iulog,"('iondrag: btesla=',   /,(6e12.4))") btesla
    ! write(iulog,"('iondrag: bdown=',    /,(6e12.4))") bdown(:,lchnk)
    ! write(iulog,"('iondrag: beast=',    /,(6e12.4))") beast(:,lchnk)
    ! write(iulog,"('iondrag: bnorth=',   /,(6e12.4))") bnorth(:,lchnk)
    ! write(iulog,"('iondrag: omega_o2p=',/,(6e12.4))") omega_o2p
    ! write(iulog,"('iondrag: omega_op=' ,/,(6e12.4))") omega_op
    ! write(iulog,"('iondrag: omega_nop=',/,(6e12.4))") omega_nop

    !-------------------------------------------------------------------------------
    ! Ion-neutral momentum transfer collision frequency coefficients:
    !-------------------------------------------------------------------------------
    do k = 1,pver
       do i = 1,ncol
          tnti(i,k) = 0.5_r8*(tI(i,k) + tn(i,k))           ! ave of tn & ti
          sqrt_tnti = sqrt( tnti(i,k) )
          wrk       = log10( tnti(i,k) )
          !-------------------------------------------------------------------------------
          ! Collision frequency coefficients with O2 (cm3/s):
          !-------------------------------------------------------------------------------
          rnu_o2p_o2(i,k) = 2.59e-11_r8*sqrt_tnti &        ! O2+ ~ O2 (resonant)
               *(1._r8 - .073_r8*wrk)**2
          rnu_op_o2 (i,k) = 6.64e-10_r8                    ! O+  ~ O2
          rnu_nop_o2(i,k) = 4.27e-10_r8                    ! NO+ ~ O2
          !-------------------------------------------------------------------------------
          ! Collision frequency coefficients with O (cm3/s):
          !-------------------------------------------------------------------------------
          rnu_o2p_o(i,k) = 2.31e-10_r8                     ! O2+ ~ O
          rnu_op_o (i,k) = 3.67e-11_r8*sqrt_tnti  &        ! O+  ~ O (resonant)
               *(1._r8 - .064_r8*wrk)**2*colfac     
          rnu_nop_o(i,k) = 2.44e-10_r8                     ! NO+ ~ O
          !-------------------------------------------------------------------------------
          ! Collision frequency coefficients with N2 (cm3/s):
          !-------------------------------------------------------------------------------
          rnu_o2p_n2(i,k) = 4.13e-10_r8                    ! O2+ ~ N2
          rnu_op_n2 (i,k) = 6.82e-10_r8                    ! O+  ~ N2
          rnu_nop_n2(i,k) = 4.34e-10_r8                    ! NO+ ~ N2
       end do
    end do

    !-------------------------------------------------------------------------------
    ! Sub set_mean_mass (mo_mean_mass.F90) returns barm(ncol,pver) in g/mole,
    !   however, set_mean_mass sometimes returns zero in top(?) four values 
    !   of the column, so barm is calculated here, see below.
    !
    ! call set_mean_mass(ncol, state%q(1,1,imozart), barm)
    !
    ! Major species and ion number densities (mmr to cm-3):
    !-------------------------------------------------------------------------------

    call set_mean_mass( ncol, state%q(:,:,imozart:), barm )

    do k = 1,pver
       do i = 1,ncol
          press(i)     = 10._r8*state%pmid(i,k) ! from Pa to dyne/cm^2
          !     barm(i,k)   = 1._r8 / (xo2(i,k)*rmi_o2 + xo1(i,k)*rmi_o1 + xn2(i,k)*rmi_n2)
          xnmbar(i,k)  = press(i)*barm(i,k)/(boltzmann*tn(i,k)) 
          o2_cm3(i,k)  = xo2(i,k)*xnmbar(i,k)*rmi_o2   ! o2 (cm-3)
          o1_cm3(i,k)  = xo1(i,k)*xnmbar(i,k)*rmi_o1   ! o  (cm-3)
          n2_cm3(i,k)  = xn2(i,k)*xnmbar(i,k)*rmi_n2   ! n2 (cm-3)
          o2p_cm3(i,k) = o2p(i,k)*xnmbar(i,k)*rmi_o2p  ! o2+ (cm-3)
          op_cm3 (i,k) = op (i,k)*xnmbar(i,k)*rmi_op   ! o+  (cm-3)
          nop_cm3(i,k) = nop(i,k)*xnmbar(i,k)*rmi_nop  ! no+ (cm-3)
!
!----------------------------------------------------------------------------------
!  Use sum of the 3 major ion number densities (as in tiegcm)
!----------------------------------------------------------------------------------
! 
          ne_sigmas(i,k)      = op_cm3(i,k) + o2p_cm3(i,k) + nop_cm3(i,k)
       end do
    end do

    !-------------------------------------------------------------------------------
    ! Multiply collision freq by neutral number density and sum for each ion:
    !
    ! rnu_o2p = [[o2p~o2]n(o2)+[o2p~o]n(o)+[o2p~n2]n(n2)]/w(o2p)
    ! rnu_op  = [[op ~o2]n(o2)+[op ~o]n(o)+[op ~n2]n(n2)]/w(op )
    ! rnu_nop = [[nop~o2]n(o2)+[nop~o]n(o)+[nop~n2]n(n2)]/w(nop)
    !-------------------------------------------------------------------------------
    do k = 1,pver
       do i = 1,ncol
          rnu_o2p(i,k) = rnu_o2p_o2(i,k)*o2_cm3(i,k) &
               + rnu_o2p_o (i,k)*o1_cm3(i,k) &
               + rnu_o2p_n2(i,k)*n2_cm3(i,k)
          rnu_op (i,k) = rnu_op_o2 (i,k)*o2_cm3(i,k) &
               + rnu_op_o  (i,k)*o1_cm3(i,k) &
               + rnu_op_n2 (i,k)*n2_cm3(i,k)
          rnu_nop(i,k) = rnu_nop_o2(i,k)*o2_cm3(i,k) &
               + rnu_nop_o (i,k)*o1_cm3(i,k) &
               + rnu_nop_n2(i,k)*n2_cm3(i,k)
          !-------------------------------------------------------------------------------
          ! Electron collision frequency (s-1):
          !-------------------------------------------------------------------------------
          sqrt_te = sqrt(tE(i,k))
          rnu_ne(i,k) = & 
               2.33e-11_r8*n2_cm3(i,k)*tE(i,k)*(1._r8 - 1.21e-4_r8*tE(i,k)) &
               + 1.82e-10_r8*o2_cm3(i,k)*sqrt_te*(1._r8 + 3.60e-2_r8*sqrt_te) &
               + 8.90e-11_r8*o1_cm3(i,k)*sqrt_te*(1._r8 + 5.70e-4_r8*tE(i,k))
       end do
    end do

    !-------------------------------------------------------------------------------
    ! Ratio of collision to gyro frequencies for o2+, o+, no+, ne:
    !-------------------------------------------------------------------------------
    do k = 1,pver
       do i = 1,ncol
          rnu_o2p(i,k) = rnu_o2p(i,k)*omega_o2p_inv(i)
          rnu_op (i,k) = rnu_op (i,k)*omega_op_inv (i)
          rnu_nop(i,k) = rnu_nop(i,k)*omega_nop_inv(i)
          rnu_ne (i,k) = rnu_ne (i,k)*omega_e_inv  (i)
       end do
    end do

    !-------------------------------------------------------------------------------
    ! Calculate pedersen and Hall conductivities (siemens/m):
    !
    ! Qe_fac: 1.e6 to convert number densities from cm-3 to m-3:
    !-------------------------------------------------------------------------------
    qe_fac(:ncol) = qe*1.e6_r8/btesla(:ncol)

    do k = 1,pver
       do i = 1,ncol
 
          !-------------------------------------------------------------------------------
          ! Pedersen conductivity (siemens/m):
          !-------------------------------------------------------------------------------
          sigma_ped(i,k) = qe_fac(i)                           &
               *((op_cm3   (i,k)*rnu_op (i,k)/(1._r8 + rnu_op (i,k)**2)) &
               + (o2p_cm3  (i,k)*rnu_o2p(i,k)/(1._r8 + rnu_o2p(i,k)**2)) &
               + (nop_cm3  (i,k)*rnu_nop(i,k)/(1._r8 + rnu_nop(i,k)**2)) &
               + (ne_sigmas(i,k)*rnu_ne (i,k)/(1._r8 + rnu_ne (i,k)**2)))
          !-------------------------------------------------------------------------------
          ! Hall conductivity (siemens/m):
          !-------------------------------------------------------------------------------
          sigma_hall(i,k) = qe_fac(i)           &
               *(ne_sigmas(i,k)/(1._r8 + rnu_ne (i,k)**2) &
               - op_cm3   (i,k)/(1._r8 + rnu_op (i,k)**2) &
               - o2p_cm3  (i,k)/(1._r8 + rnu_o2p(i,k)**2) &
               - nop_cm3  (i,k)/(1._r8 + rnu_nop(i,k)**2))
       end do
    end do

    call outfld ('ELECDEN ',ne_sigmas ,pcols,lchnk)
    call outfld ('SIGMAPED',sigma_ped ,pcols,lchnk)
    call outfld ('SIGMAHAL',sigma_hall,pcols,lchnk)

    !--------------------------------------------------------------------------------------------
    !  Save conductivities in physics buffer using pointer for access in ionosphere module
    !--------------------------------------------------------------------------------------------
    if ( waccmx_is('ionosphere') ) then 
      call pbuf_set_field(pbuf, PedConduct_idx,  sigma_ped(1:ncol,1:pver),  start=(/1,1/), kount=(/ncol,pver/) )
      call pbuf_set_field(pbuf, HallConduct_idx, sigma_hall(1:ncol,1:pver), start=(/1,1/), kount=(/ncol,pver/) )
    endif

    do k = 1,pver
       do i = 1,ncol
          wrk         = btesla(i)**2*avo*1.e-3_r8/xnmbar(i,k)
          lamda1(i,k) = sigma_ped(i,k)*wrk
          lamda2(i,k) = sigma_hall(i,k)*wrk
       end do
    end do

    call outfld ('LAMDA1',lamda1,pcols,lchnk)
    call outfld ('LAMDA2',lamda2,pcols,lchnk)

    do k = 1,pver
       do i = 1,ncol
          lxxnorot(i,k) = lamda1(i,k)
          lyynorot(i,k) = lamda1(i,k)*sin2dip(i)
          lxynorot(i,k) = lamda2(i,k)*sindip(i)
       end do
    end do

    !-------------------------------------------------------------------------------
    ! Rotate lambdas from local magnetic to geographic coordinates:
    !-------------------------------------------------------------------------------
    do k = 1,pver
       do i = 1,ncol
          lxx(i,k) = lxxnorot(i,k)*cos2dec(i) + lyynorot(i,k)*sin2dec(i)
          lyy(i,k) = lyynorot(i,k)*cos2dec(i) + lxxnorot(i,k)*sin2dec(i)
          wrk      = (lyynorot(i,k) - lxxnorot(i,k))*sindec(i)*cosdec(i)
          lyx(i,k) = lxynorot(i,k) - wrk
          lxy(i,k) = lxynorot(i,k) + wrk
       end do
    end do

    call outfld ('LXX     ',lxx,pcols,lchnk)
    call outfld ('LYY     ',lyy,pcols,lchnk)
    call outfld ('LXY     ',lxy,pcols,lchnk)
    call outfld ('LYX     ',lyx,pcols,lchnk)

    call physics_ptend_init(ptend, state%psetcols, "ion drag", lu=.true., lv=.true., ls=.true.)

    !-------------------------------------------------------------------------------
    ! Calculate ion drag tendencies and apply to neutral velocities:
    !-------------------------------------------------------------------------------
    call iondrag_tend( lchnk, ncol, state, ptend, pbuf,  &
         lxx, lyy, lxy, lyx, delt )

    !-------------------------------------------------------------------------------
    ! Calculate joule heating tendency and apply to temperature:
    !-------------------------------------------------------------------------------
    call jouleheat_tend( lchnk, ncol, state, ptend, pbuf,  &
         lxx, lyy, lxy, lyx )

  end subroutine iondrag_calc_ions

  !=========================================================================

  subroutine iondrag_calc_ghg (lchnk,ncol,state,ptend)

    use phys_grid,      only: get_rlat_all_p 
    use cam_history,    only: outfld
    use physics_types,  only: physics_ptend_init


    !
    !     This subroutine calculates ion drag using globally uniform
    !     ion drag tensor:
    !
    !                |alamxx       alamxy   | 
    !                |                      |
    !         lambda=|                      |
    !                |                      |
    !                |alamyx       alamyy   |
    !
    !     alamxx and alamxy are provided is data statements
    !     alamyy is obtaine from alamxx:
    !
    !
    !       alamyy = alamxx (sin(DIP_ANGLE))**2
    !
    !     where 
    !
    !       DIP_ANGLE = arctan(2.*tan(clat))
    !

    !--------------------Input arguments------------------------------------

    integer, intent(in) :: lchnk                   ! chunk identifier
    integer, intent(in) :: ncol                    ! number of atmospheric columns

    type(physics_state), intent(in) :: state
    type(physics_ptend ), intent(out) :: ptend
    


    !---------------------Local workspace-------------------------------------

    real(r8) :: clat(pcols)             ! latitudes(radians) for columns

    real(r8) alamyyi                                 ! ALAMYY
    real(r8) dipan                                   ! dip angle
    real(r8) dut(pcols,pver)
    real(r8) dvt(pcols,pver)

    integer i
    integer k

    !-------------------------------------------------------------------------

    if (.not.doiodrg) then
       call physics_ptend_init(ptend,state%psetcols,'none') !Initialize an empty ptend for use with physics_update
       return
    end if

    call physics_ptend_init(ptend, state%psetcols, "ion drag", lu=.true., lv=.true.)

    call get_rlat_all_p(lchnk, pcols, clat)

    !     calculate zonal wind drag
    dut(:,:)=0.0_r8
    do i=1,ncol
       do k=ntop_lev,nbot_lev
          dut(i,k)=-alamxyi(k)*state%v(i,k)-alamxxi(k)*state%u(i,k)
       enddo
    enddo

    !     calculate meridional wind drag
    dvt(:,:)=0.0_r8
    do i=1,ncol
       dipan=atan(2._r8*tan(clat(i)))
       do k=ntop_lev,nbot_lev
          alamyyi=alamxxi(k)*(sin(dipan))**2._r8
          dvt(i,k)=+alamxyi(k)*state%u(i,k)-alamyyi*state%v(i,k)
       enddo
    enddo

    do i=1,ncol
       do k=ntop_lev,nbot_lev
          ptend%u(i,k)=dut(i,k)
          ptend%v(i,k)=dvt(i,k)
       enddo
    enddo

    !  Write out tendencies
    call outfld('UIONTEND ',dut   ,pcols   ,lchnk  )
    call outfld('VIONTEND ',dvt   ,pcols   ,lchnk  )

    return
  end subroutine iondrag_calc_ghg

  !===================================================================================

  subroutine iondrag_tend( lchnk, ncol, state, ptend, pbuf,  &
       lxx, lyy, lxy, lyx, delt )

    !-------------------------------------------------------------------------------
    ! Calculate tendencies in U and V from ion drag tensors, which were
    !   calculated by sub iondrag_calc (module data lxx,lyy,lxy,lyx).
    ! This is called from sub iondrag_calc.
    !-------------------------------------------------------------------------------

    !-------------------------------------------------------------------------------
    ! dummy arguments
    !-------------------------------------------------------------------------------
    integer,intent(in)    :: lchnk                 ! current chunk index
    integer,intent(in)    :: ncol                  ! number of atmospheric columns
    real(r8), intent(in)  :: delt                  ! time step (s)
    real(r8), intent(in)  :: lxx(pcols,pver)       ! ion drag tensor
    real(r8), intent(in)  :: lyy(pcols,pver)       ! ion drag tensor
    real(r8), intent(in)  :: lxy(pcols,pver)       ! ion drag tensor
    real(r8), intent(in)  :: lyx(pcols,pver)       ! ion drag tensor
    type(physics_state), intent(in)    :: state ! Physics state variables
    type(physics_ptend), intent(inout) :: ptend ! Physics tendencies
    
    type(physics_buffer_desc), pointer :: pbuf(:)


    !-------------------------------------------------------------------------------
    ! Local variables
    !-------------------------------------------------------------------------------
    integer  :: i, k
    real(r8) :: dti
    real(r8) :: detr
    real(r8) :: us, vs
    real(r8) :: l11, l12, l21, l22
    real(r8) :: dui(pcols,pver)            ! zonal ion drag tendency
    real(r8) :: dvi(pcols,pver)            ! meridional ion drag tendency
    real(r8), pointer :: ui(:,:)           ! pointer to 3d zonal ion drift from edynamo
    real(r8), pointer :: vi(:,:)           ! pointer to 3d meridional ion drift from edynamo
    real(r8), pointer :: wi(:,:)           ! pointer to 3d vertical ion drift from edynamo

    !-------------------------------------------------------------------------------
    ! Get ion ExB drift from physics buffer (they were defined by either the exbdrift 
    !   module in chemistry (2d), or the dynamo module in dynamics dpie_coupling (3d), 
    !   depending on the switch empirical_ion_velocities. If using dynamo drifts, 
    !   they were put into pbuf by dp_coupling. If using empirical exbdrifts, then
    !   they are redundant in the vertical dimension (i.e., 2d only).
    !-------------------------------------------------------------------------------
    call pbuf_get_field(pbuf, ui_idx, ui)
    call pbuf_get_field(pbuf, vi_idx, vi)
    call pbuf_get_field(pbuf, wi_idx, wi)

    dti = 1._r8/delt
    !-------------------------------------------------------------------------------
    ! Zonal (du) and meridional (dv) wind drag, using ExB drift velocities
    ! from exbdrift module (pbuf):
    !-------------------------------------------------------------------------------
    do k = ntop_lev,nbot_lev
       do i = 1,ncol
          !-------------------------------------------------------------------------------
          ! 2/28/04 btf:
          ! Full ion-drag, using lambdas and ExB drifts. 
          !   This should succeed with bz = 0 (efield module)
          ! Runs:
          !   bz=-5, nstep=24 min, nsplit=4 (6 min dynamics): crashed after 2 days.
          !   bz=-5, nstep=24 min, nsplit=6 (4 min dynamics): 5 day run succeeded.
          ! See comments in efield module re bz < 0 (efield.F90).
          ! Ion drifts are from edynamo if use_dynamo_drifts=true, exbdrift otherwise.
          !-------------------------------------------------------------------------------
          us = ui(i,k) - state%u(i,k)
          vs = vi(i,k) - state%v(i,k)

          !-------------------------------------------------------------------------------
          ! Exclude ue,ve drift momentum source to avoid crashes when bz < 0 and 
          ! full 30 min timestep (partial ion-drag):
          !-------------------------------------------------------------------------------
          l11     = dti + lxx(i,k)
          l12     = lxy(i,k)
          l21     = -lyx(i,k)
          l22     = dti + lyy(i,k)
          detr    = dti/(l11*l22 - l12*l21)
          dui(i,k) = dti*(detr*(l12*vs - l22*us) + us)
          dvi(i,k) = dti*(detr*(l21*us - l11*vs) + vs)
       end do
    end do

    !-------------------------------------------------------------------------------
    ! Apply to model tendencies:
    !-------------------------------------------------------------------------------
    do k = ntop_lev,nbot_lev
       !-------------------------------------------------------------------------------
       ! Ion drag tendencies:
       !-------------------------------------------------------------------------------
       ptend%u(:ncol,k) = dui(:ncol,k)
       ptend%v(:ncol,k) = dvi(:ncol,k)
       !-------------------------------------------------------------------------------
       ! Turn off ion drag tendency:
       !-------------------------------------------------------------------------------
       !     ptend%u(:ncol,k) = 0._r8
       !     ptend%v(:ncol,k) = 0._r8
    end do
    do k = nbot_lev+1,pver
       dui(:ncol,k)      = 0._r8
       dvi(:ncol,k)      = 0._r8
       ptend%u(:ncol,k) = 0._r8
       ptend%v(:ncol,k) = 0._r8
    end do

!
! Ion drifts are either empirical (use_dynamo_drifts==false), or from edynamo
! (use_dynamo_drifts==true). See addfld calls in this source file.
! If empirical, the drifts will be 2d (i.e., redundant in the vertical dimension)
!
    if (empirical_ion_velocities) then
       call outfld ( 'UI', ui, pcols, lchnk )
       call outfld ( 'VI', vi, pcols, lchnk )
       call outfld ( 'WI', wi, pcols, lchnk )
    endif

    call outfld ( 'UIONTEND', dui, pcols, lchnk )      ! u ion drag tendency
    call outfld ( 'VIONTEND', dvi, pcols, lchnk )      ! v ion drag tendency

  end subroutine iondrag_tend

  !================================================================================================
  subroutine jouleheat_tend( lchnk, ncol, state, ptend, pbuf,  &
       lxx, lyy, lxy, lyx )
    !-------------------------------------------------------------------------------
    ! Calculate tendencies in T due to joule heating.
    ! This is called from sub iondrag_calc.
    !-------------------------------------------------------------------------------

    use physconst, only: pi,cpairv
    use phys_grid, only: get_rlon_p, get_rlat_p

    !-------------------------------------------------------------------------------
    ! dummy arguments
    !-------------------------------------------------------------------------------
    integer,intent(in)    :: lchnk                    ! current chunk index
    integer,intent(in)    :: ncol                     ! number of atmospheric columns
    real(r8), intent(in)  :: lxx(pcols,pver)          ! ion drag tensor
    real(r8), intent(in)  :: lyy(pcols,pver)          ! ion drag tensor
    real(r8), intent(in)  :: lxy(pcols,pver)          ! ion drag tensor
    real(r8), intent(in)  :: lyx(pcols,pver)          ! ion drag tensor
    type(physics_state), intent(in)    :: state       ! Physics state variables
    type(physics_ptend), intent(inout) :: ptend       ! Physics tendencies (inout)
    
    type(physics_buffer_desc), pointer :: pbuf(:)

    !-------------------------------------------------------------------------------
    ! Local variables
    !-------------------------------------------------------------------------------
    integer  :: k, i
    integer  :: max_ind(2)
    real(r8) :: us, vs
    real(r8) :: max_q
    real(r8) :: qjoule(pcols,pver)         ! joule heating
    real(r8) :: qout(pcols,pver)           ! temp for outfld
    real(r8), pointer :: ui(:,:)           ! pointer to pbuf
    real(r8), pointer :: vi(:,:)           ! pointer to pbuf

    logical, parameter :: debug = .false.

    !-------------------------------------------------------------------------------
    ! Get ion velocities from physics buffer (they were defined by exbdrift module)
    ! Ion velocities are 2d arrays, i.e., no vertical dimension.
    !-------------------------------------------------------------------------------
    call pbuf_get_field(pbuf, ui_idx, ui     )
    call pbuf_get_field(pbuf, vi_idx, vi     )
 
    do k = ntop_lev,nbot_lev
       !   write(iulog,"('qjoule: k=',i3,' u=',/,(6e12.4))") k,state%u(:,k)
       !   write(iulog,"('qjoule: k=',i3,' v=',/,(6e12.4))") k,state%v(:,k)
       !   write(iulog,"('qjoule: k=',i3,' lxx=',/,(6e12.4))") k,lxx(:,k)
       !   write(iulog,"('qjoule: k=',i3,' lxy=',/,(6e12.4))") k,lxy(:,k)
       !   write(iulog,"('qjoule: k=',i3,' lyx=',/,(6e12.4))") k,lyx(:,k)
       !   write(iulog,"('qjoule: k=',i3,' lyy=',/,(6e12.4))") k,lyy(:,k)
       do i = 1,ncol
          us           = ui(i,k) - state%u(i,k)
          vs           = vi(i,k) - state%v(i,k)
          qjoule(i,k)  = us*us*lxx(i,k) + us*vs*(lxy(i,k) - lyx(i,k)) + vs*vs*lyy(i,k)
          ptend%s(i,k) = qjoule(i,k)        ! joule heating tendency
       end do
       !   write(iulog,"('qjoule: k=',i3,' qjoule(:,k)=',/,(6e12.4))") k,qjoule(:,k)
    end do
    do k = nbot_lev+1,pver
       qjoule(:ncol,k)  = 0._r8
       ptend%s(:ncol,k) = 0._r8        ! no joule heating tendency
    end do

    sw_debug: if (debug) then
       max_q      = 100._r8*maxval( abs( qjoule(:ncol,ntop_lev:nbot_lev) )/state%t(:ncol,ntop_lev:nbot_lev) )
       max_ind(:) = maxloc( abs( qjoule(:ncol,ntop_lev:nbot_lev) )/state%t(:ncol,ntop_lev:nbot_lev) )
       i = max_ind(1)
       k = max_ind(2)
       if( lchnk == 25 ) then
          i = 14
          k = 3
          write(iulog,*) ' '
          write(iulog,*) '-------------------------------------------------------'
          write(iulog,*) 'jouleheat_tend: lon,lat = ',get_rlon_p(lchnk,14)*180._r8/pi, get_rlat_p(lchnk,14)*180._r8/pi
          write(iulog,*) 'jouleheat_tend: dt,t,max% dt/t = ',qjoule(i,k)/cpairv(i,k,lchnk),state%t(i,k),max_q, &
               ' @ lchnk,i,k = ',lchnk,max_ind(:)
          write(iulog,*) 'jouleheat_tend: lxx,xy,yx,yy   = ',lxx(i,k),lxy(i,k),lyx(i,k),lyy(i,k)
          write(iulog,*) 'jouleheat_tend: u,ui,v,vi      = ',state%u(i,k),ui(i,k),state%v(i,k),vi(i,k)
          write(iulog,*) 'jouleheat_tend: us,vs          = ',ui(i,k) - state%u(i,k),vi(i,k) - state%v(i,k)
          write(iulog,*) 'jouleheat_tend: du,dv          = ',ptend%u(i,k),ptend%v(i,k)
          write(iulog,*) 'jouleheat_tend: dt'
          write(iulog,'(1p,5g15.7)') qjoule(max_ind(1),ntop_lev:nbot_lev)/cpairv(max_ind(1),ntop_lev:nbot_lev,lchnk)
          write(iulog,*) '-------------------------------------------------------'
          write(iulog,*) ' '
          ! stop 'diagnostics'
       end if
    endif sw_debug

    qout(:ncol,:) = qjoule(:ncol,:)/cpairv(:ncol,:,lchnk)

    call outfld ( 'QJOULE', qout, pcols, lchnk )

  end subroutine jouleheat_tend

!==============================================================================

  subroutine iondrag_inidat(ncid_ini, pbuf2d)

    use pio,      only: file_desc_t
    use ncdio_atm,only: infld
    use infnan,   only: nan, assignment(=)
    use cam_grid_support, only : cam_grid_check, cam_grid_id, cam_grid_get_dim_names
    use physics_buffer, only : pbuf_set_field

   ! args
    type(file_desc_t), intent(inout)   :: ncid_ini    ! Initial condition file id
    type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! Physics buffer

   ! local vars
    real(r8), pointer :: ui_tmp(:,:,:)
    real(r8), pointer :: vi_tmp(:,:,:)
    real(r8), pointer :: wi_tmp(:,:,:)
    real(r8) :: nanval
    integer  :: grid_id
    character(len=4) :: dim1name, dim2name
    character(len=*), parameter :: subname='iondrag_inidat'
    logical :: found_ui, found_vi, found_wi

    allocate(ui_tmp(pcols,pver,begchunk:endchunk))
    allocate(vi_tmp(pcols,pver,begchunk:endchunk))
    allocate(wi_tmp(pcols,pver,begchunk:endchunk))

    grid_id = cam_grid_id('physgrid')
    if (.not. cam_grid_check(grid_id)) then
      call endrun(trim(subname)//': Internal error, no "physgrid" grid')
    end if
    call cam_grid_get_dim_names(grid_id, dim1name, dim2name)

    call infld( 'UI',ncid_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
         ui_tmp, found_ui, gridname='physgrid')
    call infld( 'VI',ncid_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
         vi_tmp, found_vi, gridname='physgrid')
    call infld( 'WI',ncid_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
         wi_tmp, found_wi, gridname='physgrid')

    ionvels_read_from_file = found_ui .and. found_vi .and. found_wi

    ui_idx = pbuf_get_index('UI')
    vi_idx = pbuf_get_index('VI')
    wi_idx = pbuf_get_index('WI')

    if (ionvels_read_from_file) then
       call pbuf_set_field(pbuf2d, ui_idx, ui_tmp)
       call pbuf_set_field(pbuf2d, vi_idx, vi_tmp)
       call pbuf_set_field(pbuf2d, wi_idx, wi_tmp)
    else
       nanval=nan
       call pbuf_set_field(pbuf2d, ui_idx, nanval)
       call pbuf_set_field(pbuf2d, vi_idx, nanval)
       call pbuf_set_field(pbuf2d, wi_idx, nanval)
    endif

    deallocate( ui_tmp )
    deallocate( vi_tmp )
    deallocate( wi_tmp )

    if (masterproc) then
       write(iulog,*) 'iondrag_inidat: ionvels_read_from_file = ',ionvels_read_from_file
    end if
 
  end subroutine iondrag_inidat

end module iondrag
